Integrative analysis reveals therapeutic potential of pyrviniym pamoate in Merkel cell carcinoma


Pyrvinium pamoate treated MCC cells analysis

back to home

library(biomaRt)
library(edgeR)
library(limma)
library(ggplot2)
library(ggrepel)
library(clusterProfiler)
library(org.Hs.eg.db)
library(plotly)
library(dplyr)
library(DT)
library(viper)
library(dorothea)
library(ggpubr)

Importing data

gene_count<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/read_counts.rds")
genemart<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/imr90_wgcna/genemart.rds")
gse39612_tumor.vs.normal<-read.delim2(file = "/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/gse39612/gse39612_tumor_vs_normal_DEGs_analysis.csv", sep = ",", row.names = 1)

Performing DEGs analysis on pyrvinium treated MCC cell lines

Dividing samples to WaGa and MKL1 cells separately

df_MKL1<-gene_count[,1:15]
print(head(df_MKL1))
##                   Mcon1 Mcon2 Mcon3 MD24h1 MD24h2 MD24h3 MD6h1 MD6h2 MD6h3
## ENSG00000223972.5     5     4     7      3     15      9    19     3     8
## ENSG00000227232.5   795   734   718    566    872    715   903   937   621
## ENSG00000278267.1    12     8     3      9     15      4    12     2     3
## ENSG00000243485.5     2     5     4      1      2      0     4    10     2
## ENSG00000284332.1     0     0     0      0      0      0     0     0     0
## ENSG00000237613.2     0     0     0      0      0      0     0     2     0
##                   MP24h1 MP24h2 MP24h3 MP6h1 MP6h2 MP6h3
## ENSG00000223972.5      8      7      7     4     1     4
## ENSG00000227232.5    781    818    721   618   921   770
## ENSG00000278267.1      9      7      4    11     6     2
## ENSG00000243485.5      3      1      9     1     2     0
## ENSG00000284332.1      0      0      0     0     0     0
## ENSG00000237613.2      0      0      0     0     0     0
df_WaGa<-gene_count[, 16:30]
print(head(df_WaGa))
##                   Wcon1 Wcon2 Wcon3 WD24h1 WD24h2 WD24h3 WD6h1 WD6h2 WD6h3
## ENSG00000223972.5    30    27    29     29     25     30    32    21    32
## ENSG00000227232.5   749   928  1018   1362   1215    995  1006  1274  1426
## ENSG00000278267.1     0     4     3     14     17     11     4    17    10
## ENSG00000243485.5    21     7     1     17      8      6     7    10    14
## ENSG00000284332.1     0     0     0      0      0      0     0     0     0
## ENSG00000237613.2     0     0     0      0      0      0     0     0     0
##                   WP24h1 WP24h2 WP24h3 WP6h1 WP6h2 WP6h3
## ENSG00000223972.5     22     22     14    37    41    32
## ENSG00000227232.5   1319   1690   1529  1196  1227  1040
## ENSG00000278267.1      6      7     15     5     6     9
## ENSG00000243485.5     12      8      7    22     6     8
## ENSG00000284332.1      0      0      0     0     0     0
## ENSG00000237613.2      0      0      0     0     0     1

Performing normalization and DEGs analysis for MKL1 samples

group_MKL1<-factor(c(rep("Mcon",3), rep("MD24h", 3), rep("MD6h", 3), rep("MP24h",3), rep("MP6h",3)))
ensemblid_MKL1<-rownames(df_MKL1)
ensemblid_MKL1<-unlist(lapply(ensemblid_MKL1, function(x) strsplit(x, "[.]")[[1]][1]))
gene_MKL1<- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol","chromosome_name","band","description","start_position","end_position", "entrezgene_id"),
                   values=ensemblid_MKL1,mart=genemart)

y<- DGEList(counts=df_MKL1,group=group_MKL1)
keep <- rowSums(cpm(y)>1) >= length(group_MKL1)/5 #cpm: compute counts per million or reads per kilobase per million#
y <- y[keep, ,keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
mm <- model.matrix(~0 + group_MKL1)
temp <- voom(y,mm, plot = T) ## voom + TMM
edat_MKL1 = temp$E
edat_MKL1<-as.data.frame(edat_MKL1)
fit <- lmFit(temp, mm)
efit <- eBayes(fit)
plotSA(efit, main="Final model: Mean-variance trend")

edat_MKL1[,"ensembl_gene_id"]<-unlist(lapply(rownames(edat_MKL1), function(x) strsplit(x, "[.]")[[1]][1]))
edat_MKL1_annotated<-left_join(edat_MKL1, gene_MKL1, by = "ensembl_gene_id")
edat_MKL1_annotated<-edat_MKL1_annotated[edat_MKL1_annotated$hgnc_symbol != "",]
edat_MKL1_annotated<-edat_MKL1_annotated[!duplicated(edat_MKL1_annotated$hgnc_symbol),]
edat_MKL1_annotated<-na.omit(edat_MKL1_annotated)
rownames(edat_MKL1_annotated)<-edat_MKL1_annotated$hgnc_symbol
edat_MKL1_annotated<-edat_MKL1_annotated[,c(1:15)]

design <- model.matrix(~0 + group_MKL1)
colnames(design)<-gsub("group_MKL1", "", colnames(design))
contr.matrix <- makeContrasts(
  MP6hvsMD6h = MP6h-MD6h, 
  MP24hvsMD24h = MP24h-MD24h, 
  MPvsMD = (MP6h+MP24h)-(MD6h+MD24h),
  MDvsMcon = (MD6h+MD24h)-Mcon,
  MPvsMcon= (MP6h+MP24h)-Mcon,
  levels = colnames(design))
colnames(fit$coefficients)<-gsub("group_MKL1","",colnames(fit$coefficients))
vfit <- contrasts.fit(fit, contrasts=contr.matrix)
efit <- eBayes(vfit)

DEGs_df_all_comparisons_MKL1<-list()
for (i in 1:length(colnames(efit$coefficients))){
  coef_name<-colnames(efit$coefficients)[i]
  print(coef_name)
  DEGs_df<-topTable(efit, coef = coef_name, n = Inf)
  rownames(DEGs_df)<-unlist(lapply(rownames(DEGs_df), function(x) strsplit(x, "[.]")[[1]][1]))
  DEGs_df$ensembl_gene_id<-rownames(DEGs_df)
  DEGs_df$hgnc_symbol<-gene_MKL1[match(rownames(DEGs_df), gene_MKL1$ensembl_gene_id), "hgnc_symbol"]
  DEGs_df$entrezgene_id<-gene_MKL1[match(rownames(DEGs_df), gene_MKL1$ensembl_gene_id), "entrezgene_id"]
  DEGs_df_all_comparisons_MKL1[[i]]<-DEGs_df
} #loop to generate a df list for DEGs without any filter.
names(DEGs_df_all_comparisons_MKL1)<-colnames(efit$coefficients)

Performing normalization and DEGs analysis for WaGa samples

group_Waga<-factor(c(rep("Wcon", 3), rep("WD24h", 3), rep("WD6h",3), rep("WP24h", 3), rep("WP6h",3)))
ensemblid_waga<-rownames(df_WaGa)
ensemblid_waga<-unlist(lapply(ensemblid_waga, function(x) strsplit(x, "[.]")[[1]][1]))
genes_waga <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol","chromosome_name","band","description","start_position","end_position","entrezgene_id"),
              values=ensemblid_waga,mart= genemart) 

y<- DGEList(counts=df_WaGa,group=group_Waga)
keep <- rowSums(cpm(y)>1) >= length(group_Waga)/5 #cpm: compute counts per million or reads per kilobase per million#
y <- y[keep, ,keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
mm <- model.matrix(~0 + group_Waga)
temp <- voom(y,mm, plot = T) ## voom + TMM
edat_waga = temp$E
edat_waga<-as.data.frame(edat_waga)
fit <- lmFit(temp, mm)
efit <- eBayes(fit)
plotSA(efit, main="Final model: Mean-variance trend")

edat_waga[,"ensembl_gene_id"]<-unlist(lapply(rownames(edat_waga), function(x) strsplit(x, "[.]")[[1]][1]))
edat_waga_annotated<-left_join(edat_waga, genes_waga, by = "ensembl_gene_id")
edat_waga_annotated<-edat_waga_annotated[edat_waga_annotated$hgnc_symbol != "",]
edat_waga_annotated<-edat_waga_annotated[!duplicated(edat_waga_annotated$hgnc_symbol),]
edat_waga_annotated<-na.omit(edat_waga_annotated)
rownames(edat_waga_annotated)<-edat_waga_annotated$hgnc_symbol
edat_waga_annotated<-edat_waga_annotated[,c(1:15)]

design <- model.matrix(~0 + group_Waga)
colnames(design)<-gsub("group_Waga", "", colnames(design))
contr.matrix <- makeContrasts(
  WP6hvsWD6h = WP6h-WD6h, 
  WP24hvsWD24h = WP24h-WD24h, 
  WPvsWD = (WP6h+WP24h)-(WD6h+WD24h),
  WDvsWcon = (WD6h+WD24h)-Wcon,
  WPvsWcon= (WP6h+WP24h)-Wcon,
  levels = colnames(design))
colnames(fit$coefficients)<-gsub("group_Waga","",colnames(fit$coefficients))
vfit <- contrasts.fit(fit, contrasts=contr.matrix)
efit <- eBayes(vfit)

DEGs_df_all_comparisons_waga<-list()
for (i in 1:length(colnames(efit$coefficients))){
  coef_name<-colnames(efit$coefficients)[i]
  print(coef_name)
  DEGs_df<-topTable(efit, coef = coef_name, n = Inf)
  rownames(DEGs_df)<-unlist(lapply(rownames(DEGs_df), function(x) strsplit(x, "[.]")[[1]][1]))
  DEGs_df$ensembl_gene_id<-rownames(DEGs_df)
  DEGs_df$hgnc_symbol<-genes_waga[match(rownames(DEGs_df), genes_waga$ensembl_gene_id), "hgnc_symbol"]
  DEGs_df$entrezgene_id<-genes_waga[match(rownames(DEGs_df), genes_waga$ensembl_gene_id), "entrezgene_id"]
  DEGs_df_all_comparisons_waga[[i]]<-DEGs_df
} #loop to generate a df list for DEGs without any filter.
names(DEGs_df_all_comparisons_waga)<-colnames(efit$coefficients)

Visualizing DEGs in volcano plots

DEGs_df_all_comparisons_waga<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_DEGs_all_realigned_waga.RDS")
DEGs_df_all_comparisons_MKL1<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_DEGs_all_realigned_MKL1.RDS")

generate_volcano_plot_for_Waga_cell<-function(samplename){
  volc<-data.frame(log2fc = DEGs_df_all_comparisons_waga[[samplename]][,"logFC"],
                   sig = -1*log(as.numeric(DEGs_df_all_comparisons_waga[[samplename]][,"adj.P.Val"]), base = 10),
                   genes = DEGs_df_all_comparisons_waga[[samplename]][, "hgnc_symbol"])
  volc$de <- "NS"
  volc$de[volc$sig>-1*log10(0.05)] <- "Significant"
  volc$de[abs(volc$log2fc) >= 1 & volc$sig >10]<-"label"
  volc$delabel <- volc$genes
  volc[which(volc$de != "label"), "delabel"]<- NA
  id<-samplename
  p<-ggplot(volc,aes(x=log2fc,y=sig,col=de))+
    theme_linedraw()+
    theme_light()+
    geom_point(col = ifelse(volc$de == "NS", "#bdbdbd", ifelse(volc$de =="label","#FF8000","#9ecae1")))+
    geom_text_repel(size = 2.5, colour = "black", aes(label = ifelse(delabel != "<NA>", delabel, "")), max.overlaps = 30) +
    geom_vline(xintercept=c(-1, 1), col="red",linetype=4)+
    geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
    ggtitle(as.character(id))+
    ylab("-log10(padj)") +
    xlab(paste0(id,"\n Log2(Fold Change)")) #+
  #xlim(-3,3)+
  #ylim(0,40)
  ggsave(p, width = 14, height = 10, file = paste0("/xdisk/mpadi/jiawenyang/result/merkel_cell/pp_in_mcc_paper/figure/pyrvinium_treated_", samplename, "_DEGs_volcano_plot.pdf"))
  return(p)
}

generate_volcano_plot_for_Waga_cell("WP24hvsWD24h")

generate_volcano_plot_for_MKL1_cell<-function(samplename){
  volc<-data.frame(log2fc = DEGs_df_all_comparisons_MKL1[[samplename]][,"logFC"],
                   sig = -1*log(as.numeric(DEGs_df_all_comparisons_MKL1[[samplename]][,"adj.P.Val"]), base = 10),
                   genes = DEGs_df_all_comparisons_MKL1[[samplename]][, "hgnc_symbol"])
  volc$de <- "NS"
  volc$de[volc$sig>-1*log10(0.05)] <- "Significant"
  volc$de[abs(volc$log2fc) >= 1 & volc$sig >10]<-"label"
  volc$delabel <- volc$genes
  volc[which(volc$de != "label"), "delabel"]<- NA
  id<-samplename
  p<-ggplot(volc,aes(x=log2fc,y=sig,col=de))+
    theme_linedraw()+
    theme_light()+
    geom_point(col = ifelse(volc$de == "NS", "#bdbdbd", ifelse(volc$de =="label","#FF8000","#9ecae1")))+
    geom_text_repel(size = 2.5, colour = "black", aes(label = ifelse(delabel != "<NA>", delabel, "")), max.overlaps = 30) +
    geom_vline(xintercept=c(-1, 1), col="red",linetype=4)+
    geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
    ggtitle(as.character(id))+
    ylab("-log10(padj)") +
    xlab(paste0(id,"\n Log2(Fold Change)"))
  ggsave(p, width = 14, height = 10, file= paste0("/xdisk/mpadi/jiawenyang/result/merkel_cell/pp_in_mcc_paper/figure/pyrvinium_treated_", samplename, "_DEGs_volcano_plot.pdf"))
  return(p)
}

generate_volcano_plot_for_MKL1_cell("MP24hvsMD24h")

Overlapping GSE39612 MCC tumor vs. normal skin significant DEGs with pyrvinium vs.dmso treated MCC cell lines significant DEGs and visualizing with volcano plot

WP24hvsWD24h<-DEGs_df_all_comparisons_waga$WP24hvsWD24h
WP24hvsWD24h_sig_upregulated_DEGs<-WP24hvsWD24h[WP24hvsWD24h$logFC>=1 & WP24hvsWD24h$adj.P.Val <=0.05,]
WP24hvsWD24h_sig_downregulated_DEGs<-WP24hvsWD24h[WP24hvsWD24h$logFC<=-1 & WP24hvsWD24h$adj.P.Val <=0.05,]

gse39612_tumor.vs.normal$logFC<-as.numeric(gse39612_tumor.vs.normal$logFC)
gse39612_tumor.vs.normal$adj.P.Val<-as.numeric(gse39612_tumor.vs.normal$adj.P.Val)
gse39612_sig_upregulated_DEGs<-gse39612_tumor.vs.normal[gse39612_tumor.vs.normal$logFC>=1,]
gse39612_sig_upregulated_DEGs<-gse39612_sig_upregulated_DEGs[gse39612_sig_upregulated_DEGs$adj.P.Val<=0.05,]
gse39612_sig_downregulated_DEGs<-gse39612_tumor.vs.normal[gse39612_tumor.vs.normal$logFC<=-1, ]
gse39612_sig_downregulated_DEGs<-gse39612_sig_downregulated_DEGs[gse39612_sig_downregulated_DEGs$adj.P.Val<=0.05,]

GSE39612_upregulated_DEGs_decreased_in_pyr_waga_DEGs<-intersect(rownames(gse39612_sig_upregulated_DEGs), WP24hvsWD24h_sig_downregulated_DEGs$hgnc_symbol)
GSE39612_downregulated_DEGs_increased_in_pyr_waga_DEGs<-intersect(rownames(gse39612_sig_downregulated_DEGs), WP24hvsWD24h_sig_upregulated_DEGs$hgnc_symbol)

volc<-data.frame(log2fc = DEGs_df_all_comparisons_waga[["WP24hvsWD24h"]][,"logFC"],
                 sig = -1*log10(DEGs_df_all_comparisons_waga[["WP24hvsWD24h"]][,"adj.P.Val"]),
                 genes = DEGs_df_all_comparisons_waga[["WP24hvsWD24h"]][, "hgnc_symbol"])
volc$de <- "NS"
volc$de[volc$sig>-1*log10(0.05)] <- "Significant"
annotate_down_df<-volc[volc$genes %in% GSE39612_upregulated_DEGs_decreased_in_pyr_waga_DEGs,]
annotate_up_df<-volc[volc$genes %in% GSE39612_downregulated_DEGs_increased_in_pyr_waga_DEGs,]

p<-ggplot(volc,aes(x=log2fc,y=sig))+
  theme_linedraw()+
  theme_light()+
  geom_point(col = "grey50", size = 3, alpha = 0.3)+
  #geom_text_repel(size = 2.5, colour = "black", aes(label = ifelse(delabel != "<NA>", delabel, ""))) +
  geom_point(data = annotate_down_df , # New layer containing data subset il_genes       
             size = 3,
             shape = 21,
             alpha = 0.8,
             fill = "firebrick",
             colour = "black")+
  geom_text_repel(data = annotate_down_df, aes(label = genes), size = 2.8, fontface = "bold", colour = "black", force = 1, max.overlaps = 20) +
  geom_point(data = annotate_up_df, # New layer containing data subset il_genes       
             size = 3,
             shape = 21,
             fill = "blue",
             alpha = 0.8,
             colour = "black")+
  geom_text_repel(data = annotate_up_df, aes(label = genes), size = 2.8, fontface = "bold", colour = "black", force = 1, box.padding = 0.5, max.overlaps = 20) +
  geom_vline(xintercept=c(-1, 1), col="red",linetype=4)+
  geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
  ylab("-log10(padj)") +
  xlab("WaGa pyrvinium treated vs. control \n Log2(Fold Change)")
p<-p+theme(axis.title=element_text(size=20,face="bold"),
           axis.text = element_text(size =15, face="bold"))
p

MP24hvsMD24h<-DEGs_df_all_comparisons_MKL1$MP24hvsMD24h
MP24hvsMD24h_sig_upregulated_DEGs<-MP24hvsMD24h[MP24hvsMD24h$logFC>=1 & MP24hvsMD24h$adj.P.Val <=0.05,]
MP24hvsMD24h_sig_downregulated_DEGs<-MP24hvsMD24h[MP24hvsMD24h$logFC<=-1 & MP24hvsMD24h$adj.P.Val <=0.05,]

GSE39612_upregulated_DEGs_decreased_in_pyr_MKL1_DEGs<-intersect(rownames(gse39612_sig_upregulated_DEGs), MP24hvsMD24h_sig_downregulated_DEGs$hgnc_symbol)
GSE39612_downregulated_DEGs_increased_in_pyr_MKL1_DEGs<-intersect(rownames(gse39612_sig_downregulated_DEGs), MP24hvsMD24h_sig_upregulated_DEGs$hgnc_symbol)

volc<-data.frame(log2fc = DEGs_df_all_comparisons_MKL1[["MP24hvsMD24h"]][,"logFC"],
                 sig = -1*log10(DEGs_df_all_comparisons_MKL1[["MP24hvsMD24h"]][,"adj.P.Val"]),
                 genes = DEGs_df_all_comparisons_MKL1[["MP24hvsMD24h"]][, "hgnc_symbol"])
volc$de <- "NS"
volc$de[volc$sig>-1*log10(0.05)] <- "Significant"
annotate_down_df<-volc[volc$genes %in% GSE39612_upregulated_DEGs_decreased_in_pyr_MKL1_DEGs,]
annotate_up_df<-volc[volc$genes %in% GSE39612_downregulated_DEGs_increased_in_pyr_MKL1_DEGs,]

p<-ggplot(volc,aes(x=log2fc,y=sig))+
  theme_linedraw()+
  theme_light()+
  geom_point(col = "grey50", size = 3, alpha = 0.3)+
  #geom_text_repel(size = 2.5, colour = "black", aes(label = ifelse(delabel != "<NA>", delabel, ""))) +
  geom_point(data = annotate_down_df , # New layer containing data subset il_genes       
             size = 3,
             shape = 21,
             alpha = 0.8,
             fill = "firebrick",
             colour = "black")+
  geom_text_repel(data = annotate_down_df, aes(label = genes), size = 2.8, fontface = "bold", colour = "black", force = 1, max.overlaps = 20) +
  geom_point(data = annotate_up_df, # New layer containing data subset il_genes       
             size = 3,
             shape = 21,
             fill = "blue",
             alpha = 0.8,
             colour = "black")+
  geom_text_repel(data = annotate_up_df, aes(label = genes), size = 2.8, fontface = "bold", colour = "black", force = 1, box.padding = 0.5, max.overlaps = 20) +
  geom_vline(xintercept=c(-1, 1), col="red",linetype=4)+
  geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
  ylab("-log10(padj)") +
  xlab("MKL1 pyrvinium treated vs. control \n Log2(Fold Change)")
p<-p+theme(axis.title=element_text(size=20,face="bold"),
           axis.text = element_text(size =15, face="bold"))
p

Performing GO-over representation analysis on DEGs from different cell lines and time of treatment

DEGs_df_all_comparisons_waga<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_DEGs_all_realigned_waga.RDS")
WaGa_6h_up<-DEGs_df_all_comparisons_waga$WP6hvsWD6h[DEGs_df_all_comparisons_waga$WP6hvsWD6h$adj.P.Val <= 0.05 & 
                                                      DEGs_df_all_comparisons_waga$WP6hvsWD6h$logFC >= 1,]
WaGa_6h_down<-DEGs_df_all_comparisons_waga$WP6hvsWD6h[DEGs_df_all_comparisons_waga$WP6hvsWD6h$adj.P.Val <= 0.05 & DEGs_df_all_comparisons_waga$WP6hvsWD6h$logFC <= -1,]
WaGa_6h_up[,"Cell_line"]<-rep("WaGa", nrow(WaGa_6h_up))
WaGa_6h_up[,"time"]<-rep("6h", nrow(WaGa_6h_up))
WaGa_6h_up[,"group"]<-rep("up", nrow(WaGa_6h_up))

WaGa_6h_down[,"Cell_line"]<-rep("WaGa", nrow(WaGa_6h_down))
WaGa_6h_down[,"time"]<-rep("6h", nrow(WaGa_6h_down))
WaGa_6h_down[,"group"]<-rep("down", nrow(WaGa_6h_down))

WaGa_24h_up<-DEGs_df_all_comparisons_waga$WP24hvsWD24h[DEGs_df_all_comparisons_waga$WP24hvsWD24h$adj.P.Val <= 0.05 & DEGs_df_all_comparisons_waga$WP24hvsWD24h$logFC >=1, ]
WaGa_24h_down<-DEGs_df_all_comparisons_waga$WP24hvsWD24h[DEGs_df_all_comparisons_waga$WP24hvsWD24h$adj.P.Val <= 0.05 & DEGs_df_all_comparisons_waga$WP24hvsWD24h$logFC <=-1, ]

WaGa_24h_up[,"Cell_line"]<-rep("WaGa", nrow(WaGa_24h_up))
WaGa_24h_up[,"time"]<-rep("24h", nrow(WaGa_24h_up))
WaGa_24h_up[,"group"]<-rep("up", nrow(WaGa_24h_up))

WaGa_24h_down[,"Cell_line"]<-rep("WaGa", nrow(WaGa_24h_down))
WaGa_24h_down[,"time"]<-rep("24h", nrow(WaGa_24h_down))
WaGa_24h_down[,"group"]<-rep("down", nrow(WaGa_24h_down))


DEGs_df_all_comparisons_MKL1<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_DEGs_all_realigned_MKL1.RDS")
MKL1_6h_up<-DEGs_df_all_comparisons_MKL1$MP6hvsMD6h[DEGs_df_all_comparisons_MKL1$MP6hvsMD6h$adj.P.Val<=0.05 & 
                                                      DEGs_df_all_comparisons_MKL1$MP6hvsMD6h$logFC >=1 , ]
MKL1_6h_down<-DEGs_df_all_comparisons_MKL1$MP6hvsMD6h[DEGs_df_all_comparisons_MKL1$MP6hvsMD6h$adj.P.Val<=0.05 & 
                                                        DEGs_df_all_comparisons_MKL1$MP6hvsMD6h$logFC <= -1 , ]
MKL1_6h_up[,"Cell_line"]<-rep("MKL1", nrow(MKL1_6h_up))
MKL1_6h_up[,"time"]<-rep("6h", nrow(MKL1_6h_up))
MKL1_6h_up[,"group"]<-rep("up", nrow(MKL1_6h_up))

MKL1_6h_down[,"Cell_line"]<-rep("MKL1", nrow(MKL1_6h_down))
MKL1_6h_down[,"time"]<-rep("6h", nrow(MKL1_6h_down))
MKL1_6h_down[,"group"]<-rep("down", nrow(MKL1_6h_down))

MKL1_24h_up<-DEGs_df_all_comparisons_MKL1$MP24hvsMD24h[DEGs_df_all_comparisons_MKL1$MP24hvsMD24h$adj.P.Val<=0.05 & DEGs_df_all_comparisons_MKL1$MP24hvsMD24h$logFC >=1 , ]
MKL1_24h_down<-DEGs_df_all_comparisons_MKL1$MP24hvsMD24h[DEGs_df_all_comparisons_MKL1$MP24hvsMD24h$adj.P.Val<=0.05 & DEGs_df_all_comparisons_MKL1$MP24hvsMD24h$logFC <= -1 , ]
MKL1_24h_up[,"Cell_line"]<-rep("MKL1", nrow(MKL1_24h_up))
MKL1_24h_up[,"time"]<-rep("24h", nrow(MKL1_24h_up))
MKL1_24h_up[,"group"]<-rep("up", nrow(MKL1_24h_up))

MKL1_24h_down[,"Cell_line"]<-rep("MKL1", nrow(MKL1_24h_down))
MKL1_24h_down[,"time"]<-rep("24h", nrow(MKL1_24h_down))
MKL1_24h_down[,"group"]<-rep("down", nrow(MKL1_24h_down))

pyr_treated_DEGs<-list(WaGa_6h_up, WaGa_24h_up, MKL1_6h_up, MKL1_24h_up, WaGa_6h_down, WaGa_24h_down, MKL1_6h_down, MKL1_24h_down)
pyr_treated_sig_DEGs<-Reduce(rbind, pyr_treated_DEGs)
pyr_treated_sig_DEGs<-na.omit(pyr_treated_sig_DEGs)
pyr_treated_sig_FC<-pyr_treated_sig_DEGs$logFC
names(pyr_treated_sig_FC)<-pyr_treated_sig_DEGs$entrezgene_id
pyr_treated_sig_DEGs<-pyr_treated_sig_DEGs[,c("entrezgene_id", "Cell_line", "time", "group", "logFC")]
pyr_go_result<-compareCluster(entrezgene_id~Cell_line+time,
                          data = pyr_treated_sig_DEGs,
                          fun = enrichGO,
                          OrgDb = "org.Hs.eg.db",
                          ont = "BP",
                          pAdjustMethod = "BH",
                          pvalueCutoff  = 0.05,
                          qvalueCutoff  = 0.25
                          )

Visualizing the GO Terms over-representative analysis

pyr_go_result<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/pyr_treated_go_terms_enrichment.rds")
GO_24H<-pyr_go_result@compareClusterResult[pyr_go_result@compareClusterResult$time == "24h", ]
datatable(data = pyr_go_result@compareClusterResult)
require(DOSE)
require(enrichplot)
require(viridis)

pyr_go_result@compareClusterResult$Cell_line = factor(pyr_go_result@compareClusterResult$Cell_line, levels = c("WaGa", "MKL1"))

p = dotplot(pyr_go_result, 
            x = "time", 
            font.size=22, 
            size = "Count",
            by = "p.adjust",
            title = "Top Gene Ontolgy terms of pyrvinium Treated MCC Cell Lines") + 
  facet_grid(~Cell_line) + 
  scale_color_viridis(option = "viridis", direction = -1) +
  scale_fill_viridis(option = "viridis", direction = -1) +
  scale_x_discrete(limits = c("6h", "24h"))+
  scale_size_area(max_size = 20)+
  theme(strip.text.x = element_text(size = 22, face = "bold"),
        axis.text.x = element_text(size = 22, face = "bold"),
        axis.title = element_text(size = 22, face = "bold"),
        text = element_text(size = 22, face = "bold"),
        plot.title = element_text(hjust = 0.5, size = 30))
p

Performing KEGG pathway analysis on pyrvinium treated MCC cell lines

KEGG pathway analysis of UP regulated genes by pyrvinium treatment (24 hours)

#For MKL1 24h ==================
pyr_treated_MKL1_24h_sig_DEGs_up<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "up" & pyr_treated_sig_DEGs$Cell_line == "MKL1" &   pyr_treated_sig_DEGs$time == "24h",]
pyr_treated_MKL1_24h_sig_FC_up<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "up" & pyr_treated_sig_DEGs$Cell_line == "MKL1" & 
                                                       pyr_treated_sig_DEGs$time == "24h","logFC"]
names(pyr_treated_MKL1_24h_sig_FC_up)<-pyr_treated_MKL1_24h_sig_DEGs_up$entrezgene_id

pyr_kegg_MKL1_24h_result_up<-enrichKEGG(pyr_treated_MKL1_24h_sig_DEGs_up$entrezgene_id,
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.05,
                                qvalueCutoff  = 0.25
)

#print(head(pyr_kegg_MKL1_24h_result_up@result))

edox_up_MKL1_24h<-setReadable(pyr_kegg_MKL1_24h_result_up, "org.Hs.eg.db", "ENTREZID")

p_MKL1_24h_up<-enrichplot::cnetplot(edox_up_MKL1_24h,
             foldChange = pyr_treated_MKL1_24h_sig_FC_up,
             layout = "kk",
             color_category = "#c6dbef",
             cex_category = 2,
             showCategory = 8,
             cex_label_category = 1.2,
             cex_label_gene = 1.2,
             categortSize = "pvalue")+
  scale_color_gradient2(name='fold change', midpoint = 4,  high='#de2d26', mid='#fc9272', low='#fff5f0')
p_MKL1_24h_up<- p_MKL1_24h_up + theme(strip.text.x = element_text(size = 22, face = "bold"),
                text = element_text(size = 22, face = "bold"))
pdf(file = "/xdisk/mpadi/jiawenyang/result/merkel_cell/pp_in_mcc_paper/figure/pyrvinium_treated_MKL1_24h_upregulated_kegg_cnetplot.pdf", height = 8, width = 12)
p_MKL1_24h_up
dev.off()
## png 
##   2
#For WaGa 24h ================
pyr_treated_WaGa_24h_sig_DEGs_up<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "up" & pyr_treated_sig_DEGs$Cell_line == "WaGa" &   pyr_treated_sig_DEGs$time == "24h",]
pyr_treated_WaGa_24h_sig_FC_up<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "up" & pyr_treated_sig_DEGs$Cell_line == "WaGa" & 
                                                       pyr_treated_sig_DEGs$time == "24h","logFC"]
names(pyr_treated_WaGa_24h_sig_FC_up)<-pyr_treated_WaGa_24h_sig_DEGs_up$entrezgene_id

pyr_kegg_WaGa_24h_result_up<-enrichKEGG(pyr_treated_WaGa_24h_sig_DEGs_up$entrezgene_id,
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.05,
                                qvalueCutoff  = 0.25
)

#print(head(pyr_kegg_WaGa_24h_result_up@result))

edox_up_WaGa_24h<-setReadable(pyr_kegg_WaGa_24h_result_up, "org.Hs.eg.db", "ENTREZID")

p_WaGa_24h_up<-enrichplot::cnetplot(edox_up_WaGa_24h,
             foldChange = pyr_treated_WaGa_24h_sig_FC_up,
             layout = "kk",
             color_category = "#c6dbef",
             cex_category = 2,
             showCategory = 8,
             cex_label_category = 1.2,
             cex_label_gene = 1.2,
             categortSize = "pvalue")+
  scale_color_gradient2(name='fold change', midpoint = 4,  high='#de2d26', mid='#fc9272', low='#fff5f0')
p_WaGa_24h_up<- p_WaGa_24h_up + theme(strip.text.x = element_text(size = 22, face = "bold"),
                text = element_text(size = 22, face = "bold"))

p_WaGa_24h_up

KEGG pathway analysis of DOWN regulated genes by pyrvinium treatment (24 hours)

#For MKL1 24h ==================
pyr_treated_MKL1_24h_sig_DEGs_down<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "down" & pyr_treated_sig_DEGs$Cell_line == "MKL1" &   pyr_treated_sig_DEGs$time == "24h",]
pyr_treated_MKL1_24h_sig_FC_down<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "down" & pyr_treated_sig_DEGs$Cell_line == "MKL1" & 
                                                       pyr_treated_sig_DEGs$time == "24h","logFC"]
names(pyr_treated_MKL1_24h_sig_FC_down)<-pyr_treated_MKL1_24h_sig_DEGs_down$entrezgene_id

pyr_kegg_MKL1_24h_result_down<-enrichKEGG(pyr_treated_MKL1_24h_sig_DEGs_down$entrezgene_id,
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.05,
                                qvalueCutoff  = 0.25
)

edox_down_MKL1_24h<-setReadable(pyr_kegg_MKL1_24h_result_down, "org.Hs.eg.db", "ENTREZID")

p_MKL1_24h_down<-enrichplot::cnetplot(edox_down_MKL1_24h,
             foldChange = pyr_treated_MKL1_24h_sig_FC_down,
             layout = "kk",
             color_category = "#edf8e9",
             cex_category = 2,
             showCategory = 8,
             cex_label_category = 1.2,
             cex_label_gene = 1.2,
             categortSize = "pvalue")+
  scale_color_gradient2(name='fold change', midpoint = -4, high='#deebf7', mid='#9ecae1', low='#08519c')
p_MKL1_24h_down<- p_MKL1_24h_down + theme(strip.text.x = element_text(size = 22, face = "bold"),
                text = element_text(size = 22, face = "bold"))

p_MKL1_24h_down

#For WaGa 24h ================
pyr_treated_WaGa_24h_sig_DEGs_down<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "down" & pyr_treated_sig_DEGs$Cell_line == "WaGa" &   pyr_treated_sig_DEGs$time == "24h",]
pyr_treated_WaGa_24h_sig_FC_down<-pyr_treated_sig_DEGs[pyr_treated_sig_DEGs$group == "down" & pyr_treated_sig_DEGs$Cell_line == "WaGa" & 
                                                       pyr_treated_sig_DEGs$time == "24h","logFC"]
names(pyr_treated_WaGa_24h_sig_FC_down)<-pyr_treated_WaGa_24h_sig_DEGs_down$entrezgene_id

pyr_kegg_WaGa_24h_result_down<-enrichKEGG(pyr_treated_WaGa_24h_sig_DEGs_down$entrezgene_id,
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.05,
                                qvalueCutoff  = 0.25
)

edox_down_WaGa_24h<-setReadable(pyr_kegg_WaGa_24h_result_down, "org.Hs.eg.db", "ENTREZID")

p_WaGa_24h_down<-enrichplot::cnetplot(edox_down_WaGa_24h,
             foldChange = pyr_treated_WaGa_24h_sig_FC_down,
             layout = "kk",
             color_category = "#edf8e9",
             cex_category = 2,
             showCategory = 8,
             cex_label_category = 1.2,
             cex_label_gene = 1.2,
             categortSize = "pvalue")+
  scale_color_gradient2(name='fold change', midpoint = -3,  high='#deebf7', mid='#9ecae1', low='#08519c')
p_WaGa_24h_down<- p_WaGa_24h_down + theme(strip.text.x = element_text(size = 22, face = "bold"),
                text = element_text(size = 22, face = "bold"))

p_WaGa_24h_down

Visualizing the KEGG pathway genes with cnetplot for both cell lines

pyr_kegg_result<-compareCluster(entrezgene_id~Cell_line+group,
                                data = pyr_treated_sig_DEGs,
                                fun = enrichKEGG,
                                pAdjustMethod = "BH",
                                pvalueCutoff  = 0.05,
                                qvalueCutoff  = 0.25
)

datatable(data = pyr_kegg_result@compareClusterResult)
edox<-setReadable(pyr_kegg_result, "org.Hs.eg.db", "ENTREZID")

cnt_enrichment <- enrichplot::cnetplot(edox,
                                       cex_gene = 1.5,
                                       cex_label_gene = 1.2,
                                       cex_label_category = 2,
                                       colorEdge = T,
                                       )+
  theme(legend.text=element_text(size=10, face = "bold"),
        legend.title = element_text(size=12, face = "bold"),
        text = element_text(face = "bold", color = "black"))+
  scale_fill_manual(values=c("MKL1.down"="#46B8DA", "MKL1.up"="#D43F3A", "WaGa.down"="#357EBD","WaGa.up" = "#EEA236"))+
  scale_colour_manual(values=c("MKL1.down"="#46B8DA", "MKL1.up"="#D43F3A", "WaGa.down"="#357EBD", "WaGa.up" = "#EEA236"), 
                      labels=c("MKL1.down", "MKL1.up", "WaGa.down", "WaGa.up"))

cnt_enrichment

Running viper to predict TFs activity on pyrvinium treated MCC RNAseq data and performing two samples ttest between conditions.

prior = dorothea_hs[dorothea_hs$confidence %in% c("A","B","C"),]
regulon = df2regulon(prior)
#WaGa cell treated with DMSO and pyrvinium 
expr = readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_waga_realigned_normalized_edat.RDS")
eset = expr[,-c(1:3)]
vpres = viper(eset = eset, regulon = regulon, 
              minsize = 10, nes = T, method = 'none', 
              eset.filter = T, pleiotropy = F, verbose = F)

res = rowTtest(x = vpres[,7:12],y = vpres[,1:6])
padj = p.adjust(res$p.value,method = "BH")
names(padj) = rownames(res$p.value)

WaGa_pyr_vs_DMSO_vpres<-data.frame(TFs = rownames(res$statistic),`t_statistics` = res$statistic, `p.adjust` = padj, sign = sign(res$statistic))
WaGa_pyr_vs_DMSO_vpres[WaGa_pyr_vs_DMSO_vpres$sign == "-1","sign"]<-"down"
WaGa_pyr_vs_DMSO_vpres[WaGa_pyr_vs_DMSO_vpres$sign == "1","sign"]<-"up"

df<-na.omit(WaGa_pyr_vs_DMSO_vpres)
cols<-c("blue","firebrick")

  p1 = ggplot(data = df, 
              aes(x=reorder(TFs, log10(p.adjust)), y=-log10(p.adjust), color = sign, size = abs(t_statistics)*2),
              label = gene)+
    geom_point()+
    scale_color_manual(values = cols)+
    geom_text_repel(data = df, aes(label = TFs), size = sqrt(abs(df$t_statistics))*3, fontface = "bold", colour = ifelse(df$sign == "up", 'firebrick', 'blue'), force = 3, max.overlaps = 50) +
    ylab(paste0("-log10(p.adjust)"))+
    xlab("Transcription Factors")+
    scale_x_discrete(expand = c(0.15,0))+
    lims(y = c(0, 5.5))+
    ggtitle("Transcritpion factors activity -  Pyrvinium vs. DMSO in WaGa cells ")+
    theme(axis.text.x = element_blank(),
          legend.position="none",
          plot.title = element_text(size = 30, face = "bold"),
          panel.background = element_blank(),
          axis.line.y = element_line(),
          #plot.margin = unit(c(1,1,1,1), "cm"),
          axis.title.x = element_text(size=30, face="bold"),
          axis.title.y = element_text(size=30, face="bold"),
          axis.text.y = element_text(size=30, face="bold"))
  
  
#MKL1 cell treated with DMSO and pyrvinium 
expr = readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc/Pyr_treatment_MKL1_realigned_normalized_edat.RDS")
eset = expr[,-c(1:3)]
vpres = viper(eset = eset, regulon = regulon, 
              minsize = 10, nes = T, method = 'none', 
              eset.filter = T, pleiotropy = F, verbose = F)

res = rowTtest(x = vpres[,7:12],y = vpres[,1:6])
padj = p.adjust(res$p.value,method = "BH")
names(padj) = rownames(res$p.value)

MKL1_pyr_vs_DMSO_vpres<-data.frame(TFs = rownames(res$statistic),`t_statistics` = res$statistic, `p.adjust` = padj, sign = sign(res$statistic))

MKL1_pyr_vs_DMSO_vpres[MKL1_pyr_vs_DMSO_vpres$sign == "-1","sign"]<-"down"
MKL1_pyr_vs_DMSO_vpres[MKL1_pyr_vs_DMSO_vpres$sign == "1","sign"]<-"up"

df<-na.omit(MKL1_pyr_vs_DMSO_vpres)
cols<-c("blue","firebrick")

  p2 = ggplot(data = df, 
              aes(x=reorder(TFs, log10(p.adjust)), y=-log10(p.adjust), color = sign, size = abs(t_statistics)*2),
              label = gene)+
    geom_point()+
    scale_color_manual(values = cols)+
    geom_text_repel(data = df, aes(label = TFs), size = sqrt(abs(df$t_statistics))*3, fontface = "bold", colour = ifelse(df$sign == "up", 'firebrick', 'blue'), force = 3, max.overlaps = 50) +
    ylab(paste0("-log10(p.adjust)"))+
    xlab("Transcription Factors")+
    scale_x_discrete(expand = c(0.15,0))+
    lims(y = c(0, 3))+
    ggtitle("Transcritpion factors activity -  Pyrvinium vs. DMSO in MKL1 cells ")+
    theme(axis.text.x = element_blank(),
          legend.position="none",
          plot.title = element_text(size = 30, face = "bold"),
          panel.background = element_blank(),
          axis.line.y = element_line(),
          #plot.margin = unit(c(1,1,1,1), "cm"),
          axis.title.x = element_text(size=30, face="bold"),
          axis.title.y = element_text(size=30, face="bold"),
          axis.text.y = element_text(size=30, face="bold"))

Viper_TF_activity_plot<-ggarrange(p1, p2,
              ncol = 1,
              nrow = 2
) 
Viper_TF_activity_plot

Utilizing ARACNe-AP + VIPER to detect TF-activity in MCC datasets with MCC specific network

Preprocessing data as input for ARACNe-AP

#Read GSE39612 data ============================================
edat<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/gse39612/renormalized_samples_edat_no_celllines.RDS")
sample_list_renormalized<-readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/gse39612/renormalized_samples_list_no_celllines.RDS")
GSE39612_annotation<-readRDS(file = "/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/gse39612/GSE39612_annotation.RDS")

MCC_tumor_matrix<-edat[ ,colnames(edat) %in% c(rownames(GSE39612_annotation[GSE39612_annotation$`tissue:ch1` == "Merkel cell carcinoma",]))]
colnames(MCC_tumor_matrix)<-paste0("Sample", seq(1:30))

MCC_tumor_matrix[, "gene"]<-rownames(MCC_tumor_matrix)
MCC_tumor_matrix<-MCC_tumor_matrix[, c(31, 1:30)]
#Read GSE39612 data ============================================
MCC_tumor_matrix_vp<-edat[ ,colnames(edat) %in% c(rownames(GSE39612_annotation[GSE39612_annotation$`tissue:ch1` == "Merkel cell carcinoma",]))]
MCC_tumor_matrix_vp<-MCC_tumor_matrix_vp[, colnames(MCC_tumor_matrix_vp) %in% c(rownames(GSE39612_annotation[GSE39612_annotation$`merkel cell polyomavirus status (dna and rna):ch1`=="positive",]))]

colnames(MCC_tumor_matrix_vp)<-paste0("Sample", seq(1:13))
MCC_tumor_matrix_vp[, "gene"]<-rownames(MCC_tumor_matrix_vp)
MCC_tumor_matrix_vp<-MCC_tumor_matrix_vp[, c(14, 1:13)]
#!/bin/bash

DATA_PATH="/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne"
OUT_PATH="/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output"
ARACNE_PATH="/xdisk/mpadi/jiawenyang/bin/aracne/ARACNe-AP/dist"


java -Xmx5G -jar ${ARACNE_PATH}/aracne.jar -e ${DATA_PATH}/mcc_tumor_matrix.txt  -o ${OUT_PATH} --tfs ${DATA_PATH}/tfs-hugo.txt --pvalue 1E-8 --seed 1 \
--calculateThreshold

for i in {1..100}
do
echo "performing ${i} of 100"
java -Xmx5G -jar ${ARACNE_PATH}/aracne.jar -e ${DATA_PATH}/mcc_tumor_matrix.txt  -o ${OUT_PATH} --tfs ${DATA_PATH}/tfs-hugo.txt --pvalue 1E-8 --seed $i
done

java -Xmx5G -jar ${ARACNE_PATH}/aracne.jar -o ${OUT_PATH} --consolidate
#!/bin/bash

#SBATCH --job-name=ARACNe-AP_MCC_TUMOR_ALL
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=2
#SBATCH --mem-per-cpu=16gb
#SBATCH --time=4:00:00
#SBATCH --partition=standard
#SBATCH --account=mpadi
#SBATCH --mail-type=ALL
#SBATCH --mail-user=jiawenyang@arizona.edu
date

# load required modules
module load ant

# need to change the following path to your file
cd /xdisk/mpadi/jiawenyang/bin/aracne


#Rscript for r files
sh run_aracne.sh

date
network_mod<-read.table("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output/network.txt", header = T)
network_mod<-network_mod[,c(1:3)]
write.table(network_mod, file ="/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output/network_mod.txt", sep = "\t", quote = F, col.names = F, row.names = F)

eset<-as.matrix(MCC_tumor_matrix[,c(2:31)])
regulon<-aracne2regulon("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output/network_mod.txt", eset, format = "3col")

Applying ARACNe-AP to build TF-targets network for MCC virus positive tumor samples

network_mod<-read.table("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output_vp/network.txt", header = T)
network_mod<-network_mod[,c(1:3)]
write.table(network_mod, file ="/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output_vp/network_mod.txt", sep = "\t", quote = F, col.names = F, row.names = F)

eset<-as.matrix(MCC_tumor_matrix_vp[,c(2:14)])
regulon<-aracne2regulon("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/TF_activity/aracne/output_vp/network_mod.txt", eset, format = "3col")
## number of iterations= 29

Performing VIPER on WaGa pyrvinium treated samples

DATA_DIR<-c("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/pyr_treated_mcc")
FILES<-list.files(DATA_DIR, pattern = ".RDS", recursive = T, full.names = T)

#WaGa cell treated with DMSO and pyrvinium 
expr = readRDS(paste0(DATA_DIR, "/Pyr_treatment_waga_realigned_normalized_edat.RDS"))
eset = expr[,-c(1:3)]
eset = eset[, c(1:3, 7:9)]

vpres = viper(eset = eset, regulon = regulon, 
              minsize = 10, nes = T, method = 'none', 
              eset.filter = T, pleiotropy = F, verbose = F)

#WaGa pyrvinium vs.DMSO treated group.
res = rowTtest(x = vpres[,4:6],y = vpres[,1:3])
padj = p.adjust(res$p.value,method = "BH")
names(padj) = rownames(res$p.value)

WaGa_pyr_vs_DMSO_vpres<-data.frame(TFs = rownames(res$statistic),`t_statistics` = res$statistic, `p.adjust` = padj, sign = sign(res$statistic))
WaGa_pyr_vs_DMSO_vpres[WaGa_pyr_vs_DMSO_vpres$sign == "-1","sign"]<-"down"
WaGa_pyr_vs_DMSO_vpres[WaGa_pyr_vs_DMSO_vpres$sign == "1","sign"]<-"up"

df<-na.omit(WaGa_pyr_vs_DMSO_vpres)
cols<-c("blue","firebrick")

#label gene of interests 

label_genes<-c("NEUROD1","ASCL1","POU4F3","SOX10", "SOX11","HES1","HES6","SOX2","ATOH1",
               "TCF3", "TCF4", "TCF7", "TCF7L1", "TCF7L2", "LEF1")

df_label<-df[label_genes,]
df_label<-na.omit(df_label)

p1 = ggplot(data = df, 
              aes(x=reorder(TFs, log10(p.adjust)), y=-log10(p.adjust), color = sign, size = abs(t_statistics)*2),
              label = gene)+
    geom_point()+
    scale_color_manual(values = cols)+
    geom_point()+
    geom_label_repel(data = df_label, aes(label=df_label$TFs), size = 6, box.padding = 2, label.padding = 0.25, max.overlaps=Inf)+
    ylab(paste0("-log10(p.adjust)"))+
    xlab("Transcription Factors")+
    scale_x_discrete(expand = c(0.15,0))+
    lims(y = c(0, 5.5))+
    geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
    ggtitle("Transcritpion factors activity -  Pyrvinium vs. DMSO in WaGa cells ")+
    theme(axis.text.x = element_blank(),
          legend.position="none",
          plot.title = element_text(size = 30, face = "bold"),
          panel.background = element_blank(),
          axis.line.y = element_line(),
          #plot.margin = unit(c(1,1,1,1), "cm"),
          axis.title.x = element_text(size=30, face="bold"),
          axis.title.y = element_text(size=30, face="bold"),
          axis.text.y = element_text(size=30, face="bold"))
p1

Performing VIPER on MKL1 pyrvinium treated samples

#MKL1 cell treated with DMSO and pyrvinium 
expr = readRDS(paste0(DATA_DIR, "/Pyr_treatment_MKL1_realigned_normalized_edat.RDS"))
eset = expr[,-c(1:3)]
eset = eset[, c(1:3, 7:9)]

vpres = viper(eset = eset, regulon = regulon, 
              minsize = 10, nes = T, method = 'none', 
              eset.filter = T, pleiotropy = F, verbose = F)

#MKL1 pyrvinium vs.DMSO treated group.
res = rowTtest(x = vpres[,4:6],y = vpres[,1:3])
padj = p.adjust(res$p.value,method = "BH")
names(padj) = rownames(res$p.value)

MKL1_pyr_vs_DMSO_vpres<-data.frame(TFs = rownames(res$statistic),`t_statistics` = res$statistic, `p.adjust` = padj, sign = sign(res$statistic))

MKL1_pyr_vs_DMSO_vpres[MKL1_pyr_vs_DMSO_vpres$sign == "-1","sign"]<-"down"
MKL1_pyr_vs_DMSO_vpres[MKL1_pyr_vs_DMSO_vpres$sign == "1","sign"]<-"up"

df<-na.omit(MKL1_pyr_vs_DMSO_vpres)
cols<-c("blue","firebrick")

#label gene of interests 

label_genes<-c("NEUROD1","ASCL1","POU4F3","SOX10", "SOX11","HES1","HES6","SOX2","ATOH1",
               "TCF3", "TCF4", "TCF7", "TCF7L1", "TCF7L2", "LEF1")

df_label<-df[label_genes,]
df_label<-na.omit(df_label)

p2 = ggplot(data = df, 
              aes(x=reorder(TFs, log10(p.adjust)), y=-log10(p.adjust), color = sign, size = abs(t_statistics)*2),
              label = gene)+
    geom_point()+
    scale_color_manual(values = cols)+
    geom_point()+
    geom_label_repel(data = df_label, aes(label=df_label$TFs), size = 6, box.padding = 2, label.padding = 0.25, max.overlaps=Inf)+
    ylab(paste0("-log10(p.adjust)"))+
    xlab("Transcription Factors")+
    scale_x_discrete(expand = c(0.15,0))+
    lims(y = c(0, 5.5))+
    geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
    ggtitle("Transcritpion factors activity -  Pyrvinium vs. DMSO in MKL1 cells ")+
    theme(axis.text.x = element_blank(),
          legend.position="none",
          plot.title = element_text(size = 30, face = "bold"),
          panel.background = element_blank(),
          axis.line.y = element_line(),
          #plot.margin = unit(c(1,1,1,1), "cm"),
          axis.title.x = element_text(size=30, face="bold"),
          axis.title.y = element_text(size=30, face="bold"),
          axis.text.y = element_text(size=30, face="bold"))
p2

Performing VIPER on IMR90 ER 48h and GFP 48h samples

expr = readRDS("/xdisk/mpadi/jiawenyang/data/merkel_cell/pp_in_mcc_paper/imr90_processed/IMR90_ER_GFP_normalized_edat.RDS")
eset = expr[,c(25:27,54:56)]
vpres = viper(eset = eset, regulon = regulon, 
              minsize = 10, nes = T, method = 'none', 
              eset.filter = T, pleiotropy = F, verbose = F)

#IMR90 cell with ER and GFP
res = rowTtest(x = vpres[,4:6],y = vpres[,1:3])
padj = p.adjust(res$p.value,method = "BH")
names(padj) = rownames(res$p.value)

IMR90_ER_vs_GFP_vpres<-data.frame(TFs = rownames(res$statistic),`t_statistics` = res$statistic, `p.adjust` = padj, sign = sign(res$statistic))

IMR90_ER_vs_GFP_vpres[IMR90_ER_vs_GFP_vpres$sign == "-1","sign"]<-"down"
IMR90_ER_vs_GFP_vpres[IMR90_ER_vs_GFP_vpres$sign == "1","sign"]<-"up"

df<-na.omit(IMR90_ER_vs_GFP_vpres)
cols<-c("blue","firebrick")

label_genes<-c("NEUROD1","ASCL1","POU4F3","SOX10", "SOX11","HES1","HES6","SOX2","ATOH1",
               "TCF3", "TCF4", "TCF7", "TCF7L1", "TCF7L2", "LEF1")

df_label<-df[label_genes,]
df_label<-na.omit(df_label)

p3 = ggplot(data = df, 
              aes(x=reorder(TFs, log10(p.adjust)), y=-log10(p.adjust), color = sign, size = abs(t_statistics)*2),
              label = gene)+
    geom_point()+
    scale_color_manual(values = cols)+
    geom_point()+
    geom_label_repel(data = df_label, aes(label=df_label$TFs), size = 6, box.padding = 2, label.padding = 0.25, max.overlaps=Inf)+
    ylab(paste0("-log10(p.adjust)"))+
    xlab("Transcription Factors")+
    scale_x_discrete(expand = c(0.15,0))+
    lims(y = c(0, 5.5))+
    geom_hline(yintercept=-1*log10(0.05),col="black", linetype=4)+
    ggtitle("Transcritpion factors activity -  IMR90 ER 48h samples vs. IMR90 GFP 48h samples ")+
    theme(axis.text.x = element_blank(),
          legend.position="none",
          plot.title = element_text(size = 30, face = "bold"),
          panel.background = element_blank(),
          axis.line.y = element_line(),
          #plot.margin = unit(c(1,1,1,1), "cm"),
          axis.title.x = element_text(size=30, face="bold"),
          axis.title.y = element_text(size=30, face="bold"),
          axis.text.y = element_text(size=30, face="bold"))
p3

back to top

back to home